#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>

#include "CMesh.h"

using namespace std;

 vector<CPoint> &CMesh::getPoints()
{
	return m_P;
}

 vector<CElement>&CMesh::getElements()
{
	return m_T;
}

 vector<CPoint> &CMesh::getPbTest()
{
	return m_Points;
}

vector<CElement> &CMesh::getTbTest()
{
	return m_Elements;
}
vector<CPoint> &CMesh::getPbTrial()
{
	return m_PbTrial;
}

vector<CElement> &CMesh::getTbTrial()
{
	return m_TbTrial;
}

vector<vector<boundary_edge>> &CMesh::getAllBoundEdges()
{
	return m_allBoundaryEdges;
}
vector<vector<boundary_node>> &CMesh::getAllBoundNodes()
{
	return m_allBoundaryNodes;
}

bool CMesh::readMeshFromFile(string filename)
{
	string suffixStr = filename.substr(filename.find_last_of('.') + 1);//get the suffix of the mesh file
	if (suffixStr == "neu")
	{
		int node_number;
		int ele_number;

		ifstream f(filename, ios::in);;
		if (!f.is_open())
		{
			cout << "\n error: open mesh file error\n";
			return false;
		}
		else
		{
			try
			{
				int ele_type;
				string s = "";
				for (int i = 0; i < 6; ++i)//skip 6 line for the file header
				{
					getline(f, s);
				}
				int NGRPS, NBSETS, NDFCD, NDFVL;
				f >> node_number >> ele_number >> NGRPS >> NBSETS >> NDFCD >> NDFVL;

				for (int i = 0; i < 3; ++i)//As >> operator does not handle the '\r' in the last line
				{
					getline(f, s);
				}

				m_Points.resize(node_number);
				m_Elements.resize(ele_number);

				//	double t1;
				for (int i = 0; i < node_number; ++i)
				{
					int index;
					double x, y;
					f >> index >> x >> y;//read the node coordinates
					m_Points[i] = CPoint(x, y, 0);
				}
				for (int i = 0; i < 3; ++i)
				{
					getline(f, s);
				}

				for (int i = 0; i < ele_number; ++i)
				{
					int id; // no check is done to see that they are monotonously increasing
					int ntype, ndp;
					f >> id >> ntype >> ndp;
					ele_type = ELE_TYPE(ndp);
					m_Elements[i].set_ele_type(ELE_TYPE(ndp));//set the element type, it is also the basis type
					for (int j = 0; j < ndp; ++j)
					{
						int list;
						f >> list;
						m_Elements[i].append(list - 1);
					}
				}

				getline(f, s);
				getline(f, s);


				char str1[30];
				int nTempBoundNum;
				char charTemp[10];
				char charTepm1[10];
				int GroupNum = 0;
				int temp1, temp2, temp3;
				for (int i = 0; i < NGRPS; ++i)
				{
					f >> charTemp >> str1 >> charTepm1;
					string str = str1;
					if (str == "GROUP")
					{
						f >> str1 >> temp1 >> charTepm1 >> nTempBoundNum >> charTemp >> temp2 >> charTemp >> temp3;
						for (int i = 0; i < 3; ++i)
						{
							getline(f, s);
						}
						int p[10];
						for (int i = 0; i < nTempBoundNum / 10; ++i)
						{
							f >> p[0] >> p[1] >> p[2] >> p[3] >> p[4] >> p[5] >> p[6] >> p[7] >> p[8] >> p[9];

						}
						int nTempGroupRow = nTempBoundNum - nTempBoundNum % 10;
						for (int i = 0; i < nTempBoundNum % 10; ++i)
						{
							int temp;
							f >> temp;
						}
						for (int i = 0; i < 3; ++i)
						{
							getline(f, s);
						}
						GroupNum++;
					}
				}

				int ntemp = 0;
				int data_type;//0=node,1=element
				for (int i = 0; i < NBSETS; ++i)
				{
					f >> str1 >> data_type >> nTempBoundNum >> temp2 >> temp3;


					if (data_type == 0)//the boundary is recorded by nodes
					{
						vector<boundary_node> boundarynodes;
						for (int i = 0; i < nTempBoundNum; ++i)
						{
							boundary_node bound_node;
							int index;
							f >> index;
							bound_node.glo_node_index = index - 1;
							boundarynodes.push_back(bound_node);
						}
						m_allBoundaryNodes.push_back(boundarynodes);
					}
					else if (data_type == 1)//the boundary is recorded by element
					{
						vector<boundary_edge> boundaryedges;
						for (int i = 0; i < nTempBoundNum; ++i)
						{
							boundary_edge boun_edge;
							int index;
							f >> index >> temp1 >> temp2;//read element node on the boundary
							boun_edge.ele_index = index - 1;
							switch (temp2)
							{
							case 1:
								if (ELE_3 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[0]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[1]);
								}
								else if (ELE_6 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[0]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[2]);
								}
								else if (ELE_4 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[0]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[1]);
								}
								else if (ELE_9 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[0]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[2]);
								}

								break;
							case 2:
								if (ELE_3 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[1]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[2]);
								}
								else if (ELE_6 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[2]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[4]);
								}
								else if (ELE_4 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[1]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[2]);
								}
								else if (ELE_9 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[2]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[4]);
								}
								break;
							case 3:
								if (ELE_3 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[2]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[0]);
								}
								else if (ELE_6 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[4]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[0]);
								}
								else if (ELE_4 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[2]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[3]);
								}
								else if (ELE_9 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[4]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[6]);
								}
								break;
							case 4://case 4 only exists when the element has four edges

								if (ELE_4 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[3]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[0]);
								}
								else if (ELE_9 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[6]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[0]);
								}
								break;
							default:
								break;
							}
							boundaryedges.push_back(boun_edge);
						}
						m_allBoundaryEdges.push_back(boundaryedges);
					}
					ntemp += 1;
					for (int i = 0; i < 3; ++i)
					{
						getline(f, s);
					}
				}

				f.close();
				generatePT();

				return true;
			}
			catch (exception& e)
			{
				cout << e.what() << endl;
				return false;
			}


		}
	}
	else if (suffixStr == "msh")
	{

	}
	else
	{
		cout << "\nUnknown mesh file" << endl;
		return false;
	}
	return false;

}



void CMesh::assignBoundNodesValue(vector<boundary_node>& bound, int index_vec,BOUNDARY_TYPE type, double(*value)(CPoint,double t))
{
	for (unsigned int i = 0; i < bound.size(); ++i)
	{
		bound[i].type = type;
		bound[i].value[index_vec]=(value);
		bound[i].mark.push_back(index_vec);
	}
}

void CMesh::assignBoundEdgesValue(vector<boundary_edge>& bound, int index_vec,BOUNDARY_TYPE type, double(*value)(CPoint,double t),double(*r)(CPoint,double t))
{
	if (Neumann==type)
	{
		for (unsigned int i = 0; i < bound.size(); ++i)
		{
			bound[i].type = Neumann;
			bound[i].value[index_vec]=value;//the value
			bound[i].mark.push_back(index_vec);
		}
	}
	else if (Robin==type)	
	{
		for (unsigned int i = 0; i < bound.size(); ++i)
		{
			bound[i].type = Robin;
			bound[i].value[index_vec]=(value);//the value
			bound[i].r[index_vec]=(r);//the r value
			bound[i].mark.push_back(index_vec);
		}
	}
	
}


void CMesh::generatePT()
{

	//generate P and T from Pb_test and Tb_test
	//Pb_test and Tb_test are read from gambit .neu file
	for (unsigned int i=0;i<m_Elements.size();++i)
	{
		if (ELE_3== m_Elements[i].get_ele_type())
		{
			//The P has duplicate nodes stored, but it doesn't affect the computation because we only need the right coordinates
			CElement T;
			for (int j = 0; j < m_Elements[i].get_ele_node_num(); j = j+1)
			{
				m_P.push_back(m_Points[m_Elements[i].getIds()[j]]);
				T.append(j+i*3);
			}
			T.set_ele_type(ELE_3);
			m_T.push_back(T);
		}

		else if ( ELE_4 == m_Elements[i].get_ele_type())
		{
			//The P has duplicate nodes stored, but it doesn't affect the computation because we only need the right coordinates
			CElement T;
			for (int j = 0; j < m_Elements[i].get_ele_node_num(); j = j + 1)
			{
				m_P.push_back(m_Points[m_Elements[i].getIds()[j]]);
				T.append(j + i * 4);
			}
			T.set_ele_type(ELE_4);
			m_T.push_back(T);
		}
		
		else if (ELE_6 == m_Elements[i].get_ele_type())
		{
			//The P has duplicate nodes stored, but it doesn't affect the computation because we only need the right coordinates
			CElement T;
			for (int j = 0; j < m_Elements[i].get_ele_node_num(); j = j + 2)
			{
				m_P.push_back(m_Points[m_Elements[i].getIds()[j]]);
				T.append(j / 2 + i * 3);
			}
			T.set_ele_type(ELE_3);
			m_T.push_back(T);
		}
		else if (ELE_9 == m_Elements[i].get_ele_type())
		{
			//The P has duplicate nodes stored, but it doesn't affect the computation because we only need the right coordinates
			CElement T;
			for (int j = 0; j < m_Elements[i].get_ele_node_num() - 1; j = j + 2)
			{
				m_P.push_back(m_Points[m_Elements[i].getIds()[j]]);
				T.append(j / 2 + i * 4);
			}
			T.set_ele_type(ELE_4);
			m_T.push_back(T);
		}
	}


	//Pb_trial and Tb_trial
	m_PbTrial = m_Points;
	m_TbTrial = m_Elements;
	
}